clear
%% Problem
d = 2;
x_bound = repmat([-6,6],[d,1]);
OptSet.sol = [3,2;
    -2.805118,3.131312;
    -3.779310,-3.283186;
    3.584428,-1.848126];
OptSet.fval = 0;
OptSet.radius = 0.1;
precision_init = 0.4; % 3.848
MaximalBudget = 50000;
precision = [0.074, 0.024, 0.0074, 0.0024];

epsilon = [0.1,0.01,0.001,0.0001];

%% problem
Problem.Dimension = d;
Problem.Domain = reshape(x_bound',[1,d*2]);
Problem.Sampling = @Sampling_BoxConstraint;
NewRegionNum = 2;

%% Algorithm
AlgorithmM.n0 = 5;
AlgorithmM.QuantileLevel = 0.2;
AlgorithmM.NewBudget = 3;
AlgorithmM.MaxSampleSize = 5*(d+1);
AlgorithmM.StopCriteria = [1,MaximalBudget];
% ClearRadius0 = (x_bound(:,2)-x_bound(:,1))./(NewRegionNum.^(MaxPartitionDepth_d));
% ClearRadius = sqrt(sum(ClearRadius0.^2)); % the maximum distance within a non-partitionable region
AlgorithmM.radius = 0;
%% OPTIMIZATION
rep = 100;
TotalBudget = zeros(rep,length(epsilon));
OptNoFound = zeros(rep,length(epsilon));
Times = zeros(rep,length(epsilon));
%%
for iepsilon = 1:length(epsilon)
    OptSet.epsilon = epsilon(iepsilon);
    % problem
    MaxPartitionDepth_d = ceil(log((x_bound(:,2)-x_bound(:,1))./precision(iepsilon))./log(NewRegionNum));
    MaxPartitionDepth = sum(MaxPartitionDepth_d);
    Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
    for i=1:rep
        disp(['epsilon: ',num2str(epsilon(iepsilon)),', replication: ',num2str(i),'.'])
        tic
        [ optima, ~, SampleSet, ~, OptSetRemain ] = PRS_MMO_test_v2( Problem, AlgorithmM, OptSet );
        Times(i,iepsilon) = toc;
        disp(['time used: ',num2str(Times(i,iepsilon)),'.'])
        TotalBudget(i,iepsilon) = size(SampleSet,1);
        OptNoFound(i,iepsilon) = size(OptSetRemain,1);
    end
end
%%
clear i iepsilon OptSetRemain SampleSet optima
save('Himmelblau.mat')
a=[mean(TotalBudget)',var(TotalBudget)'];
